function [Ss, rs] = find_spillovers_ss(Ss, Hs, rs, Prob, tol,max_It)
    global bet0 bet1 tau a b w e pi_theta N_a N_theta N_w N_e w_max w_min alpha myslope...
        xi xi1 ub lb pi_e H w theta Int_pi_w sigmae elasts Pr log_e_mu log_e_sigma phi sigma  rts r_bar fp...
        r_base base_rs w_mult shares w_bar tmp_shift renorm frac myslope gamma b_init bet0_init bet1_init tau_init base_rs_init share unemp scale idx_pe
   
 
    % Rescaling each step so that in the SS, mean wage is normalized to 1
    pi_w = sum(Prob, 2);
    w_bar = dot(w, pi_w);
    w_bar = scale;



    b = b * (scale / w_bar)^(1+alpha);
    bet0 = bet0*(scale /w_bar)^((1+alpha)/xi);
    bet1 = bet1*(scale /w_bar)^((1+alpha)/xi);
    tau = tau*(scale /w_bar);
    base_rs = base_rs*scale /w_bar;
    rs = rs*scale / w_bar;
    unemp = unemp*scale  / w_bar;
    w_bar = dot(w, pi_w);
    
    for it = 1:max_It
        % Spectral Algorithm to find market clearing prices

        x_0 = rs(1:end);
        x_0(end) = base_rs;
	    rng default;
        x_1 = x_0 + randn(size(rs,1),1)/100; 
        x_1(end) = base_rs;

        f_0 = find_rents_ss(x_0, Ss,  shares, Prob);
        f_1 = find_rents_ss(x_1, Ss,  shares, Prob);
       
        max_iter = 100;
        for i = 1:max_iter 
            
            s = x_1(1:end) - x_0(1:end);
            y = f_1(1:end) - f_0(1:end);

            x_1 = x_0;
            f_1 = f_0;

            lambda = dot(s,y)/dot(y,y);
            
            step_size = 1;
            for iter_sub = 1:80
                
                x_0_t = x_1 - lambda * f_1 * step_size;
                x_0_t(end) = base_rs;
                if sum(x_0_t(1:end-1) > 0 ) == (length(x_0_t)-1)
                    break
                else
                    step_size = step_size * 0.9;
                end
            end
            x_0 = x_0_t;
            
            x_0(end) = base_rs;
            
            f_0 = find_rents_ss(x_0, Ss,  shares, Prob);
           
            diff = sqrt(dot(x_0 - x_1, x_0 - x_1));
            diff_f = sum(abs(f_0));
            
            if (diff < 1e-6) & (diff_f < 1e-3)
                rs = x_0;
                break
            elseif i == max_iter
                disp('fminsearch')
                rs = [base_rs*1.05, base_rs*1.02,base_rs]';
                x_0 = rs(1:end-1); 
                
                minf = @(x) sum(abs(find_rents_ss([x; base_rs], Ss,  shares, Prob)));
                [minx, fval, exitflag] = fminsearch(minf, x_0 );
  
                rs = [minx; base_rs]; 
                if (fval > 1e-3 | isnan(fval) | (exitflag ~= 1)) & (it > 5 )
                    fval, minx
                   
                    error("Did not find root")
                end
            end                
        end
        
        rs(end) = base_rs;

        % Calculate Choice Probabilities
        [P_S, ~, wprimeS, ~] = choiceProb(Ss, rs);
 
        P_S = sum(bsxfun(@times, P_S, reshape(pi_theta, [1,1,1,2])),4);
        popMass = bsxfun(@times, P_S, Prob);
        
        wMass = squeeze(sum(bsxfun(@times, wprimeS, popMass), 1:2));
        popMass = squeeze(sum(sum(popMass, 1:2),4));
        
        % Calculate Spillovers
        Ssi = wMass ./ popMass;
        Ssi = Ssi.^xi1;
        
        % Test for convergence
        d_S = abs(Ssi - Ss);
        
        if sum(d_S) < 2 * tol
            break
        elseif it == max_It
            display(sum(d_S(:)))
            error("Did not converge in Spillover") 
        end
        Ss = Ssi;
     
    end
    
end